Classique

ACP(normalized.count[sharedBy3, ], clusteredGenes)

Class: pca dudi
Call: dudi.pca(df = log(data + 0.1), center = TRUE, scale = TRUE, scannf = FALSE, 
    nf = 4)

Total inertia: 24

Eigenvalues:
    Ax1     Ax2     Ax3     Ax4     Ax5 
17.2861  4.7060  0.6320  0.5028  0.2615 

Projected inertia (%):
    Ax1     Ax2     Ax3     Ax4     Ax5 
 72.025  19.608   2.633   2.095   1.089 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
  72.03   91.63   94.27   96.36   97.45 

(Only 5 dimensions (out of 24) are shown)

NULL

## PNLModels

groups <- str_split_fixed(colnames(data), "_", 2)[, 1]
co2 <- str_split_fixed(groups, "", 3)[, 1]
nitrate <- factor(str_split_fixed(groups, "", 3)[, 2])
nitrate <- relevel(nitrate, "N")
fer <- factor(str_split_fixed(groups, "", 3)[, 3])
fer = relevel(fer, "F")
covariates <- data.frame(row.names = colnames(data), co2, nitrate, fer)


DEGenes <- sharedBy3
# preparation des données
counts <- round(t(data[DEGenes, ]), 0)
plnData <- prepare_data(counts = counts, covariates = covariates)


PCA_models <- PLNPCA(Abundance ~ 1 + nitrate + fer + co2 + offset(log(Offset)), data = plnData, 
    ranks = 1:10)

 Initialization...

 Adjusting 10 PLN models for PCA analysis.
     Rank approximation = 1 
     Rank approximation = 2 
     Rank approximation = 3 
     Rank approximation = 4 
     Rank approximation = 5 
     Rank approximation = 6 
     Rank approximation = 7 
     Rank approximation = 8 
     Rank approximation = 9 
     Rank approximation = 10 
 Post-treatments
 DONE!

PCA_models

--------------------------------------------------------
COLLECTION OF 10 POISSON LOGNORMAL MODELS
--------------------------------------------------------
 Task: Principal Component Analysis
========================================================
 - Ranks considered: from 1 to 10
 - Best model (greater BIC): rank = 10 - R2 = 0.97
 - Best model (greater ICL): rank = 10 - R2 = 0.97

plot(PCA_models)

myPCA_ICL <- getBestModel(PCA_models, "ICL")
plot(myPCA_ICL, ind_cols = groups, var_cols = factor(clusteredGenes[rownames(data[DEGenes, 
    ])]))

gridExtra::grid.arrange(plot(myPCA_ICL, ind_cols = groups, map = "individual", plot = FALSE), 
    plot(myPCA_ICL, var_cols = factor(clusteredGenes[rownames(data[DEGenes, ])]), 
        map = "variable", plot = FALSE), ncol = 2)

sigma(myPCA_ICL) %>% corrplot(is.corr = FALSE)

Réseau avec PLN Network

net_norm <- PLN_network(data = data, DEGenes = sharedBy3)

 Initialization...
 Adjusting 30 PLN with sparse inverse covariance estimation
    Joint optimization alternating gradient descent and graphical-lasso
    sparsifying penalty = 1.673463 
    sparsifying penalty = 1.545729 
    sparsifying penalty = 1.427745 
    sparsifying penalty = 1.318766 
    sparsifying penalty = 1.218106 
    sparsifying penalty = 1.125129 
    sparsifying penalty = 1.039249 
    sparsifying penalty = 0.9599238 
    sparsifying penalty = 0.8866536 
    sparsifying penalty = 0.8189761 
    sparsifying penalty = 0.7564644 
    sparsifying penalty = 0.6987241 
    sparsifying penalty = 0.6453911 
    sparsifying penalty = 0.5961289 
    sparsifying penalty = 0.5506269 
    sparsifying penalty = 0.508598 
    sparsifying penalty = 0.4697772 
    sparsifying penalty = 0.4339195 
    sparsifying penalty = 0.4007988 
    sparsifying penalty = 0.3702062 
    sparsifying penalty = 0.3419487 
    sparsifying penalty = 0.315848 
    sparsifying penalty = 0.2917396 
    sparsifying penalty = 0.2694714 
    sparsifying penalty = 0.2489028 
    sparsifying penalty = 0.2299043 
    sparsifying penalty = 0.2123559 
    sparsifying penalty = 0.196147 
    sparsifying penalty = 0.1811752 
    sparsifying penalty = 0.1673463 
 Post-treatments
 DONE!

Stability Selection for PLNnetwork: 
subsampling: ++++++++++++++++++++Poisson Lognormal with sparse inverse covariance (penalty = 0.167)
==================================================================
 nb_param    loglik       BIC       ICL R_squared n_edges      EBIC pen_loglik
      955 -21692.38 -23209.91 -17592.43         1     431 -23852.85  -21788.51
 density
    0.05
==================================================================
* Useful fields
    $model_par, $latent, $var_par, $optim_par
    $loglik, $BIC, $ICL, $loglik_vec, $nb_param, $criteria
* Useful S3 methods
    print(), coef(), sigma(), vcov(), fitted(), predict(), standard_error()
* Additional fields for sparse network
    $EBIC, $density, $penalty 
* Additional S3 methods for network
    plot.PLNnetworkfit() 

V(net_norm)$color <- clusteredGenes[V(net_norm)]
plot.igraph(net_norm, vertex.size = 10, vertex.label.cex = 0.5, edge.width = 0.5)

length(E(net_norm))

[1] 431

data <- toVisNetworkData(net_norm)
visNetwork(nodes = data$nodes, edges = data$edges, height = "500px", width = "100%") %>% 
    visEdges(smooth = FALSE) %>% visPhysics(solver = "forceAtlas2Based", timestep = 1, 
    minVelocity = 10, maxVelocity = 10, stabilization = F)

GENEI3

Faut-il normaliser les données avant geneie? (pris en compte dans PLN)

library(GENIE3)
TF <- read.table("TFs_PlnTFDB.txt", h = T, sep = "\t")
TF$AGI <- str_split_fixed(TF$Protein.ID, "\\.", 2)[, 1]


for (thr in seq(0.05, 0.39, by = 0.05)) {
    print(thr)
    genie(normalized.count, regressors = intersect(TF$AGI, sharedBy3), targets = sharedBy3, 
        thr = thr)
}

[1] 0.05

[1] 0.1

[1] 0.15

[1] 0.2

[1] 0.25

[1] 0.3

[1] 0.35

net <- genie(normalized.count, regressors = intersect(TF$AGI, sharedBy3), targets = sharedBy3, 
    thr = 0.25)

data <- toVisNetworkData(net)
visNetwork(nodes = data$nodes, edges = data$edges, height = "500px", width = "100%") %>% 
    visEdges(smooth = FALSE) %>% visPhysics(solver = "forceAtlas2Based", timestep = 1, 
    minVelocity = 10, maxVelocity = 10, stabilization = F)

V(net)$color <- clusteredGenes[V(net)]
data <- toVisNetworkData(net)
visNetwork(nodes = data$nodes, edges = data$edges, height = "500px", width = "100%") %>% 
    visEdges(smooth = FALSE) %>% visPhysics(solver = "forceAtlas2Based", timestep = 1, 
    minVelocity = 10, maxVelocity = 10, stabilization = F)